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Abstract 

For a Coulomb system contained in a domain A, the dielectric 
susceptibility tensor \a is defined as relating the average polariza- 
tion in the system to a constant applied electric field, in the linear 
limit. According to the phenomenological laws of macroscopic elec- 
trostatics, xa depends on the specific shape of the domain A. In this 
paper we derive, using the methods of equilibrium statistical me- 
chanics in both canonical and grand-canonical ensembles, the shape 
dependence of xa and the corresponding finite-size corrections to the 
thermodynamic limit, for a class of general ^-dimensional (u > 2) 
Coulomb systems, of ellipsoidal shape, being in the conducting state. 
The microscopic derivation is based on a general principle: the total 
force acting on a system in thermal equilibrium is zero. The results 
are checked in the Debye-Hiickel limit. The paper is a generalization 
of a previous one [L. Samaj, J. Stat. Phys. 100:949 (2000)], dealing 
with the special case of a one-component plasma in two dimensions. 
In that case, the validity of the presented formalism has already been 
verified at the exactly solvable (dimensionless) coupling T = 2. 
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1 INTRODUCTION 



For systems with short-ranged pair interactions among constituents de- 
fined in a specifically shaped domain A, the thermodynamic limit of an 
intensive quantity does not depend in general on the shape of the domain 
A and on the conditions at the boundary dA given by the surrounding 
medium. This is no longer true in the case of macroscopic systems with 
long-ranged pair interactions. A typical example is the domain-shape de- 
pendence of the dielectric susceptibility tensor for conductors predicted by 
the phenomenological laws of electrostatics [TJ. The aim of this paper is to 
derive rigorously and precisely, using the methods of equilibrium statisti- 
cal mechanics in both canonical and grand-canonical ensembles, the shape 
dependence of the dielectric susceptibility and the corresponding finite-size 
corrections to the thermodynamic limit, for a class of general classical v- 
dimensional microscopic Coulomb systems being in the conducting state. 
The case v = 1 has special features and will not be discussed here. The 
paper is a generalization of the previous one referred to as I, which was 
devoted to the microscopic derivation of the dielectric susceptibility for the 
special case of a one-component plasma in two dimensions (2D). 

In dimension v, the Coulomb potential v at a spatial position r = 
(r , r 2 , . . . , r u ), induced by a unit charge at the origin 0, is the solution of 
the Poisson equation 

Av(r) = -sj(r) (1.1) 

where s v — 2n L 'I 2 /T{v /2) is the surface area of the ^-dimensional unit 
sphere. Explicitly, 

( -ln(r/r ) ifi/ = 2, 

«(>■) - | r 2-u (1-2) 

otherwise 

^ v — 2 

Here, r — |r| and tq is an arbitrary length scale. The corresponding force 
F(r) = — Vu(r) reads 

F « = £ (1.3) 

In a (/-dimensional space, the definition of the Coulomb potential (|1.1|) 
implies in the Fourier space the characteristic small- fc behavior v(k) oc 
1/fc 2 . This maintains many generic properties (like screening) of "real" 3D 
Coulomb systems. 

We consider general Coulomb systems consisting of M mobile species 
a = 1,...,M with the corresponding charges q a , embedded in a fixed 
uniform background of charge density p^. The most studied models are the 
one-component plasma (OCP) and the symmetric two-component plasma 
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(TCP). The OCP corresponds to M = 1 with qi = q and pb of opposite sign; 
it may be convenient to define a "background density" rib by pb = —qnb- 
The symmetric TCP corresponds to M = 2 with q\ = q,q2 = —q and 
Pb = 0. The system is contained in a domain A of specified shape with a 
smooth boundary <9A. The surrounding medium is for simplicity a vacuum 
producing no image forces. The fixed background produces the one-particle 
potential pb4>b(r) where 

4> b (r)= /dVw(lr-r'l) (1.4) 

J A 

The corresponding electric field is pf,Ef,(r) where 

E fe (r)=-V0 6 (r) = j A dV^-£j_ (1.5) 

The energy of a configuration {r^ , q ai } of the charged particles plus the 
background is 

E = ^q ai q a] v(\Yi - Tj\) + Pb^2qati<Pb(ri) + #6-6 (1-6) 

i<j i 

Since the backgroud-background interaction energy term Eb-b does not de- 
pend on the particle coordinates, its particular value is irrelevant in the cal- 
culation of particle distribution functions. In the case of point particles, for 
many-component systems with at least two oppositely charged species, the 
singularity of v(r) (|1.2|) at the origin prevents the thermodynamic stability 
against the collapse of positive-negative pairs of charges: in two dimen- 
sions for small enough temperatures, in three and higher dimensions for 
any temperature. However, in those cases, one can introduce short-range 
repulsive interactions which prevent the collapse. The derivations which 
follow allow for such interactions. 

The Coulomb system in the domain A at inverse temperature f3 will be 
considered in both canonical (fixed particle numbers) and grand canonical 
(fixed species chemical potentials) ensembles. The thermal average will be 
denoted by (• • •). In terms of the microscopic density of particles of species 
a, n a (r) — J2i S a , ai S(r — r^), the microscopic densities of the total particle 
number and charge are defined respectively by 

n(r) = ^n Q (r), p(r) = q a n a (r) (1.7) 

a ol 

At one-particle level, the total particle number and charge densities are 
given respectively by 

n(r) = (n(r)>, p(r) = (p(r)) (1.8) 
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At two-particle level, one introduces the two-body densities 

= {h a (r)h a ,(r'))-{n a (r))5 a>a ,5(v-r') (1.9) 

The corresponding Ursell functions are defined by 

C/ aa /(r,r') = n^,(r,r')-n Q (rK-(r') (1.10) 

and the truncated charge-charge structure function by 

S(r,r') = (p(r)p(r')) T 

ee (p(r)p(r'))-(p(r))(p(r')) (1.11) 

The small-A: behavior of the Fourier transform of the Coulomb potential 
gives rise to exact moment constraints for the charge structure function 
S (see review J3J). In the bulk, liniA^R.- 5a (r, r') = S(\r — r'|) obeys 
the Stillinger-Lovett screening rules 01 E] which imply the zeroth-moment 
(clcctroneutrality) condition 

J d u rS{r) = (1.12) 
and the second-moment condition 

13 [ d v r\r\ 2 S(r) = -— (1.13) 

J Sv 

For finite systems, the analog of the zeroth-moment sum rule 

/ & v rS{r,r')= [ d"r'S(r,r') = (1.14) 
Ja J a 

holds only in the canonical ensemble where it reflects the trivial fact that 
the total charge in the domain A is fixed. In the grand canonical ensem- 
ble, the system is expected to exhibit charge fluctuations 6 , in which 
case (|1.14|l does not hold. The information analogous to the bulk second- 
moment condition l|1.13(l is contained in the dielectric susceptibility tensor 
XA- Let us use the notation 

d"rrV(r) i = l,...,v (1.15) 

for the ith component of the total polarization in the system. The tensor 
XA is defined as relating the average polarization to a constant applied field 
Eo, in the linear limit: 

/>>E - = Ex^'^' (1.16) 
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The linear response theory expresses the XA-components as 





|A| 



/ d"n /d v r 2 r^S , (ri,r 2 ) (1.17) 

J A J A 

where (• ■ ■) is an average defined for Eo = 0. In the canonical ensemble 
where the sum rule l|1.14f> applies, the tensor components Xa are expressible 
in another equivalent way 

Xl=-^J A d"n J d"r 2 {r\ - ri) a S(n, r 2 ) (1.18) 

As A — > R" one might naively expect that only the diagonal components 
X % = limA-fR" Xa (* = 1> ■ ■ ■ > ^) survive and, according to the bulk second- 
moment sum rule (|1.13fl , that they tend to the uniform "Stillinger-Lovett" 
(SL) value 

(1.19) 



X i S L = -~Jd»r(rfS(r) = ±- 



which does not depend on the shape of A. This is indeed true for a 
boundary-free domain like the surface of a sphere. As is explained below, 
relation l|1.19fl no longer holds in a geometry with a boundary. 

According to phenomenological electrostatics, based on plausible but 
not rigorously justified arguments, the dielectric susceptibility x of a macro- 
scopic system is related to its dielectric constant e. For the considered 
Coulomb plasma in a conducting state, the equality e _1 = implies 

Xl = f (?I T (1-20) 

where Ta is the size-invariant but shape-dependent depolarization tensor 
with position-independent components 

Tl = ~ I d»rpP, (1.21) 
A SuJa dr'dr^ v ; 

Without any loss of generality one can choose a coordinate system in which 
Ta is diagonal, T A — T A Sij, and consequently xa is also diagonal, Xa ~ 
X\Sij . Then, Eq. (|1.2()jl takes the form 

X\ = AfT (i- 22 ) 
s v T A 

For ^-dimensional ellipsoidal domains PP which will be of interest in this 
work, the TA-components <|1.21H are expressible in an alternative form as 
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with <f>b (r) defined by 1)1. 4|) . where r is an arbitrary point in A. With regard 
to the Poisson equation (|l.lf> . the diagonal elements of Ta are constrained 
by = 1- I n the special isotropic case of ^-dimensional spheres, 

T A = \ jv and the consequent x\ — v I s v is v times Xsl 01 Eq. I|1.19|l . 

The discrepancy between the naive prediction of statistical mechanics 
i|1.19fl and phenomenological electrostatics (|1.22|) was explained in a nice se- 
ries of papers by Choquard et al. The point is that the susceptibility 
is made up of a bulk contribution, which saturates quickly to the SL value 
l|1.19fl . and of a surface contribution. The surface contribution does not 
vanish in the thermodynamic limit due to the inverse-power-law behavior of 
the charge structure function at large distances along the boundary. Sum- 
ming up both contributions one gets instead of (|1.19fl the shape-dependent 
result of macroscopic electrostatics (|1.22|l . This fact was verified on the 
2D disk geometry, in the high-temperature Debye-Huckcl limit and at the 
exactly solvable coupling T — [3q 2 = 2 of the OCP. A progress towards 
the microscopic verification of formula (|1.22|) was made in paper I. There, 
the mapping of the 2D OCP, when T is an even positive integer, onto a 
discrete ID anticommuting-field theory was used for generating a sum 
rule for the charge structure function. This sum rule comes from a specific 
unitary transformation of anticommuting variables keeping a "composite" 
form of the fermionic action. For A an elliptic domain, the sum rule con- 
firms microscopically the asymptotic formula (|1.22|) and gives a finite-size 
correction term to x\ explicitly in terms of boundary contributions. 

The underlying sum rule derived for the 2D OCP seemed to be closely 
related to the logarithmic nature of the 2D Coulomb potential. We show in 
this paper that actually the sum rule is nothing but a direct consequence 
of a general principle: the total force acting on a system in thermal equilib- 
rium is zero. Using this principle, the sum rule is generalized to an arbitrary 
^-dimensional Coulomb plasma. As the result, for A a ^-dimensional el- 
lipsoidal domain, the asymptotic formula (|1.22|) for x\ 1S reproduced and 
its leading finite-size correction is obtained, in both canonical and grand 
canonical ensembles. 

The paper is organized as follows. Section 2 is devoted to the deriva- 
tion of the crucial sum rule for an arbitrary z^-dimensional Coulomb plasma. 
Based on this sum rule, the splitting of the susceptibility into its macro- 
scopic part (|1.22(l plus a corresponding finite-size correction term is shown 
for ^-dimensional domains of ellipsoidal shape in Section 3. Section 4 
presents an analysis of the finite-size correction term, dependent on the 
particular ensemble. The formalism is documented in Section 5 on the 
Debye-Hiickel limit. The check on the exactly solvable 2D OCP at cou- 
pling r = 2 has already been done in the previous paper I. Concluding 
remarks are given in Section 6. 
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2 SUM RULES 



One of us has derived several sum rules for the 2D OCP in paper I, using 
a mapping on a fermionic field theory. Actually, these sum rules are much 
more general. In the present section, the generalization of some of these 
sum rules is obtained by simple arguments about the balance of forces or 
torques. 

Writing that the total force acting on the particles is zero, at equilib- 
rium, results into a sum rule relating their density n(r) and their charge 
density p(r): 



where the first term in the l.h.s. is the force exerted by the background, 
and the second term is the force exerted by the walls. dS = d^n where n 
is the unit vector normal to the surface element dS and directed towards 
the exterior of A. This is the generalization of Eqs. (56) of paper I. For 
simplicity, we have assumed that the particle-wall interaction is a hard one, 
such that the center of each particle feels a hard wall on dA. 

A similar sum rule is obtained by assuming that a particle of species 
ai is fixed at point r l5 and writing that the total force acting on the other 
particles vanishes. Now, the force that the fixed particle exerts on the other 
ones must also be included in the force balance, which reads 



(3 Pb / d'raEifaJ^g^nW^ra.r!)- / dS 2 £ ng> ai (r 2 , n) 



where we have used that the density of particles of species a 2 at r 2 , knowing 

(2) 

that there is a particle of species oti at ri, is n a ^ ai (r 2 , ri)/n ai (ri). If there 
are short-range interactions, they must be added to the definition 1)1.3(1 of 
the Coulomb force F. Another form of Eq. 1(2.2(1 can be obtained by using 
the first BGY equation which can be written as 

Vn ai (ri) = /3p 6 E & (ri)g Ql n Ql (ri) 



With regard to the equality F(ri — r 2 ) = — F(r 2 — ri), using Eq. 1(2.3(1 for 
the last term in the l.h.s of Eq. ((2.2(1 gives 




(2.1) 
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-/ dS 2 ^ni 2 2 ) Qi (r 2 ,r 1 )-Vn Ql (r 1 ) = (2.4) 

Finally, we multiply Eq. (|2.4II by q ai and sum on ot\, we mutiply Eq. (|2.1|) 
by Pp(r{), and we substract from each other the two resulting equations, 
with the result 

/3 Pb / d w r 2 E 6 (r 2 )S(r 2 ,r 1 )=Vp(r 1 ) + / dS 2 \ q ai U a2 

Ja JdA anoa 

(2-5) 

This is the crucial sum rule which is the generalization of Eq. (60) of paper 
I. 

Although we shall not need them in the following, let us mention that 
another class of sum rules can be obtained from the balance of torques. 
For instance, in three dimensions, writing that the total torque acting on 
the particles (due to both the background and the walls) vanishes at equi- 
librium gives the sum rule 



(3p b [ d 3 r[r x E 6 (r)]p(r) + / [dS x r]n(r) = 

J A JdA 



(2.6) 



This is the generalization of Eq. (41b) of paper I. If one particle is assumed 
to be fixed at some point, one obtains the torque analog of Eq. (|2.5() 

(3p b [ d 3 r 2 [r 2 x E 6 (r 2 )]S(r 2 |ri) = [n x V]p(n) - (2.7) 

J A 

/ [dS 2 x r 2 ] V" g Ql ^ Q2Ql (r 2 ,ri) 

This is the generalization of Eq. (45b) of paper I. 

The sum rules (41a) and (45a) of paper I can also be generalized, fol- 
lowing a method developed in refs. and \T2\. However, these general- 
izations will not be described here. 



3 DERIVATION OF THE SUSCEPTIBIL- 
ITY 

Let A be a j/-dimensional ellipsoid in the reference frame defined by the 
axes of the ellipsoid, 

i— 1 x ' 

In this reference frame both tensors xa and Ta are diagonal. For the do- 
main shape under consideration, the depolarization tensor Ta is expressible 
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as (ll.23[) and independent of the point re A, while <fi b (r) is invariant under 
the transformations r l — ► — r\ This implies that 

V 

(f> b (r) = const - y 5Z T A ( r ') 2 ( 3 - 2 ) 
i=i 

The corresponding E&(r) = — V0(r) reads 

E 6 (r) = S „^T A rV (3.3) 

i=l 

where e 1 is the unit vector along the ith axis. The components of T\ for a 
2D ellipse read 

R 2 R 1 

T ' K = Tl = WTw (3 ' 4) 

The components of Ta are more complicated functions of R 1 , R 2 , R 3 for a 
3D ellipsoid P^. In the isotropic case R 1 = R of a ^-dimensional sphere, 
T A = l/i/. 

Inserting (|3.3() into the sum rule (|2.5() . and defining the components 
dS^ = dS2 • e l , one gets for each component the equality 

PpbSvTl [ d v r 2 r 2 S(r 2 ,r 1 ) = -^ T p(r 1 ) + [ dSi, V q^U^^n) 
Ja dr i J dA aua2 

(3.5) 

We multiply both sides of l|3.5|l by r\, then integrate J A d v r\ and use the 
definition (|1.17(l of the dielectric susceptibility, for obtaining 



PbS 



,T aXa |A| = f d»rr i ^p(r) + f d - nr \ f d5 * V q ai U aaai (r 2 , ri ) 

(3.6) 

Simple algebra gives 

/ A dW |^ (r) = l A d v r-^[r i p(r)]-J^r[p(r)+p b -p b ] 

dS l r l p(r)- (Q)+p b \A\ (3.7) 



where 



d^r[p(r) +/9fc ] (3.8) 



is the microscopic total charge (including the fixed background charge) in 
the domain. Provided that p b ^ 0, Eq. i|3.6[l can be thus rewritten in the 
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final form 



Xa 



|A| 



1 

|Aj 



d v nr\ / d^(p(n)n(r 2 )) r 



OA 



(3.9) 



This is the desired splitting of the susceptibility onto its macroscopic part 
plus a finite-size correction term. 
The formula H3.9[) can be further simplified in the isotropic case of a 
^-dimensional spherical domain A with a radius R l — R and a volume 
|A| = s v R v jv. Since now the components Xa do not depend on i, we 
can consider their common value xa — Yli=i Xa/ 17 - For the ^-dimensional 
sphere it holds 



dS* 



R 



dS 2 , ^24dSi = r lC os9dS2 



(3.10) 



where 9 is the angle between and r 2 . Since (p(r 1 )n(r 2 )) T depends on 
the orientations of ri and r 2 only through their angle 9, we can choose r 2 
along the 1-axis and replace f g >dS2 by s^R 1 "^ 1 . Eq. I|3.9|) takes the form 



v 

XA = — 



(3.11) 



where R = (R, 0, . . . , 0). It is sometimes convenient to express r 1 in the 
integral on the r.h.s. of H3.11[l as r 1 = R — (R — r 1 ) and in this way to 
obtain an alternative "boundary" form of Eq. H3.11JI . 



XA = — 



<|-«3*(R)) T + i 



d ,y r(i?-r 1 )(p(r)n(R)) r 



(3-12) 

The above formalism applies to the case pb ^ 0, with no restriction 
on the use of canonical or grand-canonical ensembles. When — > 0, for 
the sake of simplicity we shall restrict ourselves to the symmetric TCP in 
a //-dimensional sphere and to only microscopic states such that the total 
charge of the system is equal to zero, Q = 0. This is either the case of 
the canonical ensemble with imposed charge neutrality, or the case of a 
restricted grand-canonical ensemble when the fixed background of charge 
— Nq is first neutralized by N opposite charges +q and then ±q charges 
are added to the system in a variable number of neutral pairs ^21- Under 
these conditions, relation (|3.12|) reduces to 



XA 



v 1 

PbSv R 



[ d"r (R - r 1 ) (p(r)ri(R)) 

J A 



(3.13) 
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where the notation (• • ■) Pb is used to emphasize that the average is taken in 
presence of the background. The background-charge density pb couples to 
particle coordinates in the Boltzmann factor exp[— f3pb J. d l/ r'0b(r')/5(r')], 
where 4>b is given in (|3.2|) . In the limit pb — > 0, the thermal average (■ ■ •) Pb 
of a microscopic quantity can be expanded around pb = 0, denoted simply 
as (•■■), using the linear response theory: 

(•••>« = <->-#% / dV4(r')(-p(r')) T + 0(^) (3.14) 

Since, due to the + «-> — charge symmetry of the TCP, (/5(r)) = and 
(p(r)n(R)) = at any point r e A, relation (|3 . 1 3|) can be rewritten in the 
Pb — ► limit as follows 




x [<p(r)p(r')n(R)) - (p(r)p(r')) <n(R))] (3.15) 

We see that for the TCP with no background, three-body densities enter 
the finite-size contribution. 

4 NON-EQUIVALENCE OF ENSEMBLES 

Although the macroscopic result for the z/-dimensional sphere xa ~ v/s v is 
the same in both the canonical and grand-canonical ensembles, the finite- 
size correction term in (|3.12|l is ensemble-dependent. 

4.1 Canonical Ensemble 

In the canonical ensemble, the microscopic total charge is fixed, Q = Q. 
Let us analyze term by term the finite-size corrections appearing on the 
r.h.s. of Eq. fTTfy 

If there is some excess charge in the domain A, due to the electrostatic 
repulsion it has tendency to move to the domain boundary dA and to create 
there a macroscopic surface charge density a = Q/\dA\. We note that, as 
a consequence, (Q)/|A| = isa/R, and it is reasonable to assume that a is 
finite. The other thermal averages in (|3.12|) are assumed to be taken for a 
fixed er. 

Since the microscopic total charge does not fluctuate, (Qn(Il)) T = 0. 

One has to be cautious when identifying the R — > oo limit of the dipole 
moment in the last term with its flat hard- wall counterpart: owing to a slow 
power-law decay of the correlations along a plain hard wall the 
limit cannot be freely interchanged with the integration. In particular, let 
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us consider in ^-dimensions a semi-infinite Coulomb plasma which occupies 
the half-space x > 0; we denote by y the set of remaining (V— 1) coordinates 
normal to x. The plane at x = is charged with the uniform surface charge 
density a. It is shown in Appendix that in dimensions v = 2, 3 the R — > oo 
limit of the considered sphere dipole moment is related to the corresponding 
flat dipole moment as follows: 



lim J d v r(R-r 1 )(p(r)n(R)) T = 2 dx x J 



dy(p(x,y)n(0,0)f 



(4.1) 

The factor 2 in this equation was first observed in paper I for the case of the 
2D OCP at coupling T — 2. Its temperature-independence is also checked 
in the Debye-Hiickel limit (see the next section). 

We conclude that in v = 2, 3 dimensions the formula for the dielectric 
susceptibility tensor of the Coulomb conductor, evaluated in the canonical 
ensemble up to the leading 1/R finite-size correction term, reads 



v v 1 
XA p 



va + 2 J dxx J dy(/5(x,y)n(0,0)) T 



(4.2) 



This result can be readily extended to the pb — > limit of the symmetric 
TCP with Q = (and, consequently, a = 0), discussed at the end of the 
previous section. Using for the truncated correlation in (|4.2|) the linear 
response (|3.14|) . now in the half-space geometry with 0&(r') = — s t/ (x') 2 /2, 
one arrives at 

XA ~ ~T~ V ^\\ dxxjdyj dx {x' f J dy' 

[(p(x,y)p(x\y')n(0,O)} - (p(x, y)p(x', y')}(n(0, 0))] (4.3) 

Although the finite-size analysis was made for the ^ = 2,3 spherical ge- 
ometries, it can be simply generalized via Eq. I|3.9I) to an arbitrary v- 
dimensional ellipsoid: the leading correction term is still of the order of 1 
over the characteristic length of the domain. 



4.2 Grand Canonical Ensemble 

The grand-canonical analysis of the finite-size corrections in H3.12J1 funda- 
mentally depends on the dimension. 

Two Dimensions. In the grand canonical ensemble, necessarily the total 
charge Q vanishes and does not fluctuate 1201 (except in a very special 
case not discussed here). This is because bringing a charged particle into 
the system from a reservoir at infinity, with a hole left in the reservoir, 
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would cost an infinite energy, and this cannot be achieved with finite fu- 
gacities. Thus, the terms (Q) and (Q«-(R)) T vanish in (|3 . 1 1 1> and l|3.12[l . 
Furthermore, l|4.1|l and l|4.3[l are still valid. 

Three Dimensions. In the grand canonical ensemble, for a finite system, 
(Q) is determined by the fugacities and does not vanish, except for special 
adjustments of these fugacities. However, (Q) is at most of order R. Indeed, 
when the sphere already carries a charge Q, the work required from bringing 
one more particle of charge q into the system from the reservoir has an 
electric part qQ/R. Therefore, with finite chemical potentials, qQ/R has 
to be finite. 

The total charge does fluctuate, with a variance such that f3(Q 2 } T = R 
in the large- R limit, and the term (Qn(R)) T in l|3.12[l does not vanish, and 
is of order 1/R as shown below. 

Indeed, considering for simplicity the case of the OCP in a 3D sphere 
of radius R, (Qn(R)) T is proportional to the total charge on the sphere 
when one of the particles of charge q is fixed on the surface. Macroscopic 
electrostatics says that, when a point charge q is at distance r > R from 
the center of a grounded sphere, it induces on it a surface charge q' = 
— (R/r)q. Thus, the total charge q+q' vanishes if r = R. However, actually, 
the "surface" charge has some microscopic thickness A of the order of the 
charge correlation length, and it is better to describe approximately the 
configuration of a particle fixed on the surface as a particle at distance R 
from the center of a sphere of radius R — A. Thus q 1 = — [(R — X)/R]q, the 
total charge q + q' is of order qX/R, and (Qn(R)) T is expected to be of 
order pX/R. 

Finally, in H3.11j) and H3.12|) . in the large- R limit, the term (Q)/|A| is 
at most of order 1/R 2 an can be discarded. But the term (Qrj(R)) T gives 
to (|3.12() a contribution of order 1/R, like the dipole integral, and both 
should be kept in the leading finite-size correction. As to (|4.1() and J4.3L 
they are still valid. 

5 DEBYE-HUCKEL THEORY 

The formulas (|3.11(l and (|4.1() will now be tested in the weak-coupling limit, 
which is described by the Debye-Huckel theory, for the general system of 
M species of point particles plus a background, in two or three dimensions. 

5.1 General Formalism 

A consistent way of deriving the Debye-Hiickel theory for a finite system 
is to start with the renormalized Mayer diagrammatic expansion (which is 
reviewed, for instance, in refs. |18j and |19p. in the grand canonical ensem- 
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ble, and to make a topological reduction, replacing the fugacities by the 
densities. The weak-coupling limit for the correlation functions is obtained 
by resumming the chain diagrams with the densities taken as constants 
n a (taking into account their position-dependence near the boundary dA 
would give corrections of higher order). This is equivalent to writing the 
Ornstein-Zernicke equations with the direct correlation functions replaced 
by —J3 times the corresponding interaction potential: 



h aia2 (r 1 ,r 2 ) = -Pq a - L q a2 v(\r 1 -r 2 \) (5.1) 

d u r 3 [-Pq ai q a3 v(\ri - r 3 |)]n Q3 /i Q 3 Q2 (r 3 , r 2 ) 



E 

C*3 



where the correlation functions h are related to the Ursell functions by 
U ai a 2 ( r i; r 2) = n ai n a2 h aia2 (ri,r2). The set (15.111 of M 2 coupled equa- 
tions can be transformed into one equation. Indeed, let us make the ansatz 
that the solution is of the form 

h aia2 (r 1 ,r 2 ) = ~f3q ai q a2 G(r 1: r 2 ) (5.2) 

Using l|5.2l) in l|5.1[) one does check that these Ornstein-Zernicke equations 
are satisfied provided that G obeys the integral equation 

G(r l! r 2 )=t;(|r 1 -r 2 |)-— f d v rsv(\ ri - r 3 |)G(r 3 , r a ) (5.3) 

s v J A 

where k 2 = s„/?^ Q n a < laj tne Debye length is 1/k, Using (|5.2(l one finds 
(p( ri )n(r 2 )) T = -^©(n.ra) + p*(n - r 2 ) (5.4) 



1 / K 2 X 



and 

5(r 1 ,r 2 )^(p(r 1 )p(r 2 )) T --^^-j G(r x , r 2 ) + j-8{v x - r 2 ) (5.5) 

The integral equation (|5.3|) for G can be transformed into a differential 
equation by taking the Laplacian with respect to ri . One obtains the usual 
Debye-Hiickel equation for the screened Coulomb potential G 

[Ax - k 2 ]G{y u r 2 ) = -sj^ - r 2 ) (5.6) 

However, in a finite system, the differential equation (|5.t)|) must be supple- 
mented by boundary conditions. In the present approach, these boundary 
conditions are provided by the integral equation (|5.3() . 

It has been seen in Section 3 that, in general, when there is a back- 
ground, the finite-size correction to the susceptibility can be expressed in 
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terms of the two-body correlation appearing in (|3.11|) . while, in the limit 
of no background, one obtains the more complicated expression l|3.15ll in 
terms of a three-body correlation. The Dcbye-Hlickel theory has the very 
special feature that this complication does not arise. Indeed, since p = —pb, 
one sees in (|5.4ll that (p(ri)fi(r2)) T / Pb is expressed in terms of the two- 
body function G even in the limit pb — > 0. Furthermore (Q) = 0. Therefore 
<|3.11[1 still involves only a two-body correlation in this limit pb — > 0. 

5.2 2D Disk 

In an infinite plane, l|5.6Jl gives G(ri,r2) = Kq[k\yi — r2 1), where Ko is a 
modified Bessel function. In a finite disk of radius R, the solution is of the 
form 

OO 

G(n,r 2 ) = J2[It(s<)Kz( S> ) + ath{sx)h{s2)]wco$l6 (5.7) 

where si t2 — ^^1,2, s< and s> are the smallest and the largest, respectively, 
of si and s 2 , Ig and Kg are modified Bessel functions, and ag a coefficient 
to be determined; pg is the Neumann factor /l*o = 1) fJ>£ = 2 for t > 1. In 
the square bracket of l|5.7l) the first term corresponds to an expansion of 
Kq(k\y\ — r2 1), while the second term corresponds to the general symmetric 
solution of H5.6(l without the r.h.s. 8 term. 

The determination of ao from the integral equation (|5.3(l has been dis- 
cussed in ref. |2D], where it has been argued that the length scale ro in the 
2D Coulomb potential v must be made infinite at the end of the calculation. 
The result is a = K x {Z)/h{Z), where Z = kR. 

For determining ag when £ > 1, we consider the integral equation l|5.3[l . 
and use for G the expansion (|5.7|l and for v the expansion 

v(\n - r 2 |) = - In |ri ~ 1-21 = - In ^ + V - ( T -A <x»l(9 2 -0i) (5.8) 

In the angular integral on #3, only the terms involving the same I in the 
two expansions (|5.7Jl and (|5.8|l survive. In terms of the square bracket in 
lfB~7jl . i.e. 

Gg{r u r 2 ) = Ig( S< )Kg( S> ) + agIg( Sl )Ig(s 2 ) (5.9) 
one obtains, when n > r 2 , 

2Gg(r 1 ,r 2 ) = \ (gf £ ^\ (^)Vs, r 2 ) 
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One could solve (|5.1U|) . However, it is simpler to remark that it implies 
dGe(ri,r 2 )/dri\ ri= ii — — (£/R)Gg(R 7 r 2 ). Therefore, using the definition 
l|5.9|l for Ge(R,r 2 ) gives 

K' l (Z)+a i I' i (Z) = --[K e (Z) + a e I e (Z)], l>\ (5.11) 

a relation that Choquard et al.jS] have obtained by another method, in- 
volving a continuation of H5.6(l outside the disk; that method led to some 
ambiguity for determining a§. From (|5.11(l . using simple relations obeyed 
by the Bessel functions, one obtains 

This final equation turns out to be valid for all £, including I = 0. 

Using H5.7JI and (|5.12|1 in Ij5.2|) . one can easily check the perfect screening 
expected in two dimensions, even in the grand canonical ensemble: the 
charge in the cloud around a particle of charge q a2 is —q a2 , 

I & 2 riy^ q ai n ai h aia2 (Yi,Y 2 ) = ~q a2 (5.13) 

We now turn to the dielectric susceptibility. From its definition (jl . 1T|> . 
using the present 5(ri,r 2 ), one obtains [5] 



1 

XA = - 
7T 



1 2h{nR) 



kRI {kR) 



(5.14) 



Alternatively, one can use the general method of the present paper and 
check the expression (|3.11|) . Here (Q) — 0, and pb — —p. Let — £>disk be 
the dipole moment defined as the integral in (|3.11() (we call this dipole mo- 
ment —D rather than D for using the same notation as in the Appendix). 
Only the part I = 1 of G contributes to this integral. Using (/5(r)n(R)) T 
from l|5.4|) and G± from (|5.9|) with ai from l|5.11|l gives, after simple ma- 
nipulations on the Bessel functions, 

2ph( K R) 

Ddisk = -^rT (5 ' 15) 

It should be remarked that, since in 2D (Qn(R)) T = 0, -Ddisk is also the 
integral in l|3.12|l . Using (|5.15() in (|3.11|) or l|3.12(l . one retrieves the same 
Xa as in Ij5.14|l . In the large- R limit, in Ij5.14|l I\ (kR) / Iq(kR) — > 1 and one 
sees that the correction term is indeed of order 1/R: 



1G 



The dipole moment D^ at for a flat wall, in the 2D Debye-Huckel theory, 
has been computed in |16|. It can be checked that, in the limit R — > oo, 
Ddisk does have twice the value found for D^t- 



5.3 3D Sphere 

In infinite space, l|5.6|l gives G(ri,r 2 ) = exp(— re|ri — r 2 |)/|ri — r 2 |. In a 
finite sphere of radius R, the same considerations as in 2D now give |5] 



G(rx,r 2 )=£ 



21 + 1 r 



^+i(s<)^ + i(s>) +hl l+ ^{si)l i+ x{s 2 ) Pe(cosl 



(5.17) 

where Pi is a Legendre polynomial. As in 2D, the coefficients bi can be 
determined by using the integral equation i|5.3[l , with the same result as in 
ref. 0: 

b -lZt^T (5 ' 18) 

With our method, there is no special problem or ambiguity with the case 
£ = 0. 

As expected, there is no perfect screening, since the starting point was 
the grand canonical ensemble. Using 1)5. 17|) and (|5.18|) in l|5.2|l gives 



d 3 ri ^<7 ai n Ql /i QlQ2 (ri,r 2 ) = -q a 



1 



sinh Kr 2 
Kr 2 cosh kR 



(5.19) 



rather than —q a2 - 

The dielectric susceptibility, computed from its definition i|1.17[l is found 
to be 



1 - 



3J| (kR) 
kRIi (kR) 



3 

i?-»oo An 



l-± 
nR 



(5.20) 



Alternatively, one can use the general method of the present paper. Again 
(Q) = 0, pb = — /O, and only the part I = 1 of G contributes to the integral 
in (|3.11|l . One retrieves for the susceptibility the result (|5.20() . 

It should be noted that, in 3D, Q fluctuates and (<5".(R)) T 7^ 0. One 
finds 

Ii (kR) n 

(Qh(R)) T =P * 7 m ~ (5.21) 
kR1_i(kR) r^oo kR 

in agreement with the qualitative estimate of Section 4.2. Therefore, with 
£> sp h defined as the integral in l|3.12[) . the equivalence of (|3.11() and (|3.12() 

gives 

D sph ~ (5.22) 
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and the finite-size correction to \A can be decomposed as 

3 

R->oo 4n 

where the term 1 / kR is the contribution from (Qn(R)) T and the term 2/kR 
is the contribution from the dipolc moment Z? sp h seen fom the boundary. 
Again, in the limit R — > oo, -D S ph docs have twice the value found for a flat 
wall in the Debye-Hiickel theory. 



kR kR 



(5.23) 



6 CONCLUSION 

Macroscopic electrostatics predicts a shape-dependent value for the dielec- 
tric susceptibility of a conductor (the response, sometimes called polariz- 
ability, to a uniform applied electric field). In the present paper, it has been 
shown that classical (i.e. non-quantum) equilibrium statistical mechanics 
of a large class of microscopic models results into a dielectric susceptibility 
which is the sum of the macroscopic value plus an explicit finite-size correc- 
tion. Thus, the limits of validity of macroscopic electrostatics are clearly 
exhibited. 

The basis for the microscopic derivation only is that the total force 
acting on a system vanishes at equilibrium. It is quite surprising that such 
a simple statement is sufficient, and the reason for that still is an open 
problem. 

Our approach deals with models of Coulomb systems made of charged 
particles embedded in a uniformly charged background. The case of no 
background is dealt with as a limiting case. It seems that our method 
cannot be used for directly starting with a system without a background. 

Classical statistical mechanics has been used. It gives an acceptable 
phenomenological description of some systems such as electrolytes or molten 
salts. We have not attempted to deal with a more fundamental description 
of real matter based on quantum statistical mechanics of point charges. 

APPENDIX: DIPOLE MOMENTS 

We briefly summarize known facts about the large-distance behavior of 
particle correlations along a plain, possibly homogeneously charged, hard- 
wall in v = 2, 3 dimensions. Let us first review the case of a semi- infinite 
Coulomb plasma which occupies the half-space x > 0; y denotes the set of 
[y — 1) coordinates normal to the x-axis. According to ref. (TS|, for the 
charge-density correlator one expects an asymptotic power-law behavior 
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along the boundary of type 

(^yHx',0)) T ~*M, | y | ... x (A . n 



where gv(x, x'), which as a function of x and x' has a fast decay away from 
the wall, obeys the relation 



dx g„(x, x') = / dxx dy(p(x, y)n(x' , 0)) T , ^ = 2,3 



•v 
2 

(A.2) 

valid for any x' > 0. A behavior of type IjA.lll at large distances was 
observed also in the large- R limit of the ^-dimensional sphere HO EI • F° r 
two points r and r' inside the sphere, it is only necessary to identify x 
and x' with the corresponding point distances from the sphere surface and 
|y| with the Euclidean distance (chord) of the point projections onto the 
sphere surface: 

x = R-r,x' = R-r'; |y| = 2i?sin(0/2) (A.3) 

where 9 is the angle between r and r'. At small distances, an infinites- 
imal deformation of a flat boundary towards the sphere has a negligible 
effect on the correlations. We can therefore write, on both microscopic and 
macroscopic scales, that, as the radius of the sphere R — ► oo, 

(pW-( 1 - , )) T | sphcrc -(^,y)n(x',o)) T | flat (A.4) 

In the dipole integral on the r.h.s. of (|3.12|) . the correlator of interest 
is taken at the point r' = R fixed at the boundary, which corresponds to 
x' = in (|A.3|I . To simplify the notation, we define 



^ sph (x,0) = (p(r)n(R)) T | sphcrc , 4> Rat (x, |y|) = (p(x, y)n(0, 0)) T | flat 

(A.5) 

Within the identification 1A.3|) with x' = 0, the asymptotic R — > oo equiv- 
alence (IA.4|I now takes the form 

^ S ph(x, 9) ~ # at (x, |y|); |y| - 2Rsm(9/2) (A.6) 

Our task is to relate the R — > oo limit of the sphere dipole moment -D S ph 
seen from the boundary and the flat dipole moment £)fl a t, defined as follows 

D sph = [ d-V(# - r^sphOR - (A.7) 

J A 

£>aat = J dxx J dyi/j &at (x, |y|) (A.8) 

Because of slight differences, the derivations of the relation are made sep- 
arately for 2D (with notation "disk" instead of "sph" ) and 3D. 
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2D Disk 

Using the substitution x = R — r and writing r 1 = (R — x)[l — 2 sin 2 (8/2)], 
the disk dipole moment (|A.7|) is expressible as 

rR I* 7T 

Aiisk = / dxx(R-x) d6ip disk (x,0) 



pR PIT 

-2 / da; {R - xf I dd si 

JO J-TT 



srr/(0/2)V>disk(M) (A.9) 

In the large- R limit, we make use of the transformation l|A.6(l to get £>disk = 
h + h , where 



1\ 
RJ 



x \ f dv 
dxx(l--) I y ^fl a t(a;,y) (A.10) 

2fi Ji_j£ 



4R 2 



'-2R I _ VI 



Since ipfiat(x,y) as a function of x has a fast decay away from the wall, 
the x/R terms in Ii and I2 can be neglected in comparison with the unity 
when R — > 00. After simple algebra, one finds for I\ 

lim L = D flat + / dxx lim 2R [ dt ( , 1 - 1^ i/> flat (a;.2iR) 



(A.12) 

Considering ip^ &t {x ,2tR) ~ g 2 (x,0)/ (2tR) 2 implies a converging integral 
over t, so that liniR^oo ii = -Dflat- As concerns the second integral I2, it 
can be analogously written as 

r°° r 1 <u 

lim I 2 = / dx lim (2i?) 2 / t 2 ^ at { Xl 2tR) (A.13) 

H-kxj Jo fl^oo < /_ 1 ^/l - £2 

As above, we consider the leading asymptotic behavior of ^>a a t(x, 2tR), 
with the result 

r°° r 1 At 

lim/ 2 =/ dx 52 (x,0) / = Dflat (A.14) 

R ^°° Jo J-i VI - i 

Here, relation (| A. 2|> with s 2 = 27r was applied at x' — 0. We conclude that 
lim D disk = 2D flat , 1/ = 2 (A.15) 

it— >00 
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3D Sphere 

In 3D, the volume element d 3 r = r 2 drdf2, where the angular part dfi = 
sin 6 d6 dtp with 8 £ (0, tt) and tp £ (0, 2n). Using the substitution x = R — r 7 
the sphere dipole moment (|A.7|I reads 

D sph = J dxx(R- x) 2 J dQip sph (x,9) 

+2 J dx(R- x) 3 J dQ sm 2 {9/2)4> sph (x,e) (A.16) 



In the large- R limit, the transformation (|A.6|) implies Z? sp h = Ix+ 12, where 
h = I dxx[l--J (2tt) ^ dyy^ &at (x,y) (A.17) 



h = 2 J dxR(l--) (2tt) ^ dyy^^ at (x,y) (A.18) 



Here, we have used J dQ — {2-k/R 2 ) J^ R dyy. As in 2D, the x/R terms are 
neglected as R — > 00. Thus, liiUR_-c Ii = £Wt and 

lim I 2 = 2tt dx lim (2R) 3 / d^ 3 ^ flat (x, 2ti?) (A.19) 

R^oo J R^oo J Q 

The leading asymptotic behavior ipR at (x, 2tR) ~ 53 (x, 0)/(2ii?) 3 as R —> 00 
gives 

/* 00 

lim J 2 = 2tt / dx 53 (x, 0) = £> flat (A.20) 
where the relation (|A.2|) with S3 = 47r was applied at x 1 = 0. Finally, 

lim A P h = 2^ flat , t/ = 3 (A.21) 

R—>oo 
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